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I.  INTRODUCTION 


The  purpose  of  this  report  is  to  present  the  results  of  an  examination  of 
the  thermodynamic  performance  of  the  discrete  injection  transpiration-cooled 
nosetip  (TCNT)  for  reentry  vehicles.  By  analyzing  the  physical  model  of  the 
system  numerically,  predictions  for  nosetip  surface  temperatures  during 
flight  may  be  obtained.  Test  data  indicate  that  in  the  case  of  a reentry  vehicle 
clear  air  trajectory,  aerodynamic  heating  alone  can  cause  the  temperature  of 
the  nosetip  to  increase  beyond  the  material's  melting  point.  In  an  environ- 
ment where  particles  of  snow  or  dust  impinge  upon  the  nosetip  surface,  it  has 
been  found  that  erosion  of  the  surface  is  more  prevalent.  The  model  which 
describes  the  system  must  therefore  take  into  account  the  possibility  of  a 
receding  surface. 

An  underlying  assumption  in  the  physical  description  of  the  problem  is 
that  the  plate-like  elements  which  compose  the  nosetip  may  be  represented  as 
individual  webs  (fins),  while  the  coolant  flowing  through  the  nosetip  flow 
passageways  may  be  considered  to  absorb  energy  from  the  individual  fins, 
with  the  coolant  then  removing  heat  by  flowing  out  of  the  system.  The  energy 
balance  described  in  Appendix  A results  in  two  coupled,  one -dimensional 
transient  heat  conduction  equations  which  involve  temperatures  of  the  fin  and 
the  fluid.  Surface  heat  fluxes,  which  become  part  of  the  fin  boundary  condi- 
tions, are  taken  to  be  those  at  the  stagnation  point  of  the  nosetip.  The  fin  is 
represented  as  a semi- infinite  slab,  thus  implying  that  the  distance  from  the 
nosetip  surface  to  the  orifice  through  which  coolant  enters  is  much  greater 
than  the  thermal  penetration  depth. 

An  implicit  finite  difference  scheme  is  used  to  solve  the  equations 
numerically  on  a computer.  Iterative  techniques  are  incorporated  to  ensure 
that  dramatic  parameter  changes  (either  in  boundary  conditions  or  in  the 
coupling  of  the  equations)  do  not  drive  the  system  unstable.  Temperature 
distributions  for  clear  air  and  weather  environments  are  compared  with  actual 
test  data  recorded  in  a ballistic  range. 
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II.  FORMULATION  OF  THE  PHYSICAL  MODEL 


The  nature  of  this  problem  is  such  that  fin  and  fluid  temperatures 
influence  one  another.  If  the  fin  temperature  Tp  is  a function  of  time  due  to 
changing  surface  conditions,  it  may  be  concluded  that  the  "bulk"  tempera- 
ture of  the  fluid  Tg  also  changes  with  time.  Separate  energy  balances  are 
performed  in  Appendix  A to  yield  two  partial  differential  equations.  For  the 
purposes  of  this  analysis,  the  fin  material  is  assumed  to  be  stainless  steel, 
and  the  fluid  is  considered  to  be  water. 


A.  ONE-DIMENSIONAL  TRANSIENT  CONDUCTION 
EQUATIONS 

The  governing  equations  for  fin  and  fluid  temperatures  as  functions  of 
time  and  distance  are  derived  in  Appendix  A and  shown  below. 

1.  FIN  EQUATION: 


aiF  a i dTF 

pfcpf  = ^ \kf  ^r. 


FLUID  EQUATION: 


;)w  + f • Q 


T PFS/  5x 


dT  d ( dTR\  / v 

/?BCPB  T)V~  ~ yKB  ~C TxTj  + (/?BCPBS  + mCPB) 


--I'Q 

x b 


whe  re 


Q = h(Tg  - Tjt-),  the  internal  convective  heat  transfer  rate  per 
unit  length 


s = the  recession  or  erosion  rate  of  fin  surface 

m = the  mass  flow  rate  per  unit  area  of  the  coolant 

The  dimensions  and  other  terms  refer  to  the  representation  of  the  physical 
model  shown  in  Fig.  1.  The  dimensions  a1  and  b'  are  substituted  for  a and  b, 
respectively,  in  the  equations  when  x lies  between  t - s dt  and  Also, 
for  this  range  of  x,  m'  should  replace  m,  where  m'  =°  brh/b1. 
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Fig.  1.  Fin  Geometry 

It  is  apparent  that  in  the  fluid,  energy  is  transported  in  the  flow  direc- 
tion by  both  convection  and  conduction.  The  Peclet  number  (DpUCp/K) 
provides  a measure  of  the  relative  importance  of  convection  to  conduction 
and,  for  the  conditions  of  interest,  has  a minimum  value  on  the  order  of 
several  hundred.  Calculations  done  to  assess  the  importance  of  axial  con- 
duction have  indicated^  that  the  error  involved  in  dropping  the  conduction 
term  is  negligible  for  values  of  the  Peclet  number  greater  than  one  hundred. 
Thus,  with  the  deletion  of  the  conduction  term,  Eq.  (2)  becomes 

aTR  / 2 

^BCPB  = (^BCPB6  + mCPB ) “c^T  ' b ’ Q (3 


Ralph  P.  Stein,  "Liquid  Metal  Heat  Transfer,  " Advances  in  Heat  Transfer, 
Vol.  3,  Academic  Press  (1966). 


This  revised  fluid  equation  simplifies  the  numerical  analysis  and  aids  in  the 
spatial  stabilization  of  the  fluid  temperatures. 


The  source  term  Q is  dependent  upon  the  type  of  heat  transfer  mecha- 
nism that  couples  the  fin  to  the  coolant.  In  this  case,  either  heat  transfer 
without  phase  change  of  the  fluid  or  subcooled  boiling  heat  transfer  is  postu- 
lated. When  no  phase  change  occurs,  that  is,  when  the  temperature  of  the 
fin  is  less  than  the  saturation  temperature  of  the  fluid,  Q takes  the  form 


Nu  . K. 


Q - — 


(tb  - tf) 


where  Nu  is  the  Nusselt  number,  and  has  been  assigned  an  average  value  of  4 
for  the  ballistic  range  and  flight  cases  reported  herein;  and  is  the  hydraulic 
diameter,  determined  as 


H 1 + b/w 


where  w is  the  width  of  the  fluid  and  fin  elements.  As  an  estimate  for  the 

subcooled  boiling  heat  transfer  rate,  an  empirical  relation  postulated  by 

. , 2 . . 
horn  is  used: 


Q(Btu/ft2-sec)  = 0.  0536  eP(T. 


T ) 

sat’ 


(temperatures  are  in  R) 


where  T^^^,  is  the  saturation  temperature  of  the  coolant  (found  as  a function 

of  the  external  pressure  and  P is  defined  to  be  P/630,  where  P is  the  external 

pressure  in  psia.  For  the  case  where  T is  greater  than  T„.  Q should  be 

r oAT 

calculated  using  both  of  the  foregoing  relations.  The  larger  value  is  then 
chosen  to  represent  the  source  term,  since  the  heat  transfer  mechanism  for 
subcooled  boiling  is  greatly  improved  over  that  without  phase  change. 

2 ~ 

J.R.S.  Thom,  et  al.,  "Boiling  in  Subcooled  Water  in  Tubes  and  Anulli,  " 
Proc.  Institute  Mech.  Eng.,  3C180:226  (1965-1966). 


It  has  been  suspected  from  prior  usat  > that  the  equation  for  subcooled 
boiling  heat  transfer  is  overly  optimistic,  i.e.,  the  equation  overpredicts 
the  amount  of  heat  removed  from  the  fin.  For  the  present,  however,  thia 
relation  shall  be  used  until  new  correlations  can  be  formulated. 

The  erosion  rate  of  the  fin  s depends  upon  the  free- stream  conditions 
outside  the  shock  wave.  An  empirical  formula  used  to  obtain  the  erosion  rate 
i s a s follows  : 


s 


1 


awe 

CN  PF 


oo 


whe  re 

P\  = the  liquid  water  content  of  the  free  stream  (the  mass  of 

water  per  volume  of  air 

U = the  free-stream  velocity,  which  is  equivalent  to  the 

velocity  of  the  vehicle 


c 


N 


an  erosion  coefficient,  herein  taken  to  be  a function  of 
the  fin  surface  temperature 


PF 


the  density  of  the  fin  material 


For  the  purposes  of  this  analysis,  it  is  assumed  that  in  a clear  air  environ- 
ment the  erosion  term  is  always  zero.  At  present,  the  formulation  does  not 
model  tin  mate ria  1 the rma Uy  removed  when  surface  temperatures  exceed  the 
material  melt  temperature.  The  important  information  to  be  obtained  in  the 
clear  air  case  is  whether  or  not  the  fin  surface  temperature  actually  reaches 
the  melt  temperature.  1 he  value  for  mass  flux  of  the  coolant  m is  assumed 
to  be  supplied  as  a function  of  time  or  altitude  of  the  vehicle.  Coolant  mass 
flux  will  be  specificed  by  design  to  restrict  surface  temperature  to  below 
material  melt  values. 


1$.  SURFACE  BOUNDARY  CONDITION  FOR  THE 
FIN  EQUATION 


The  second-order  partial  differential  equation  which  governs  the 
temperature  in  the  fin  indicates  that  two  boundary  conditions  are  needed 
for  complete  solution.  The  surface  condition  in  the  fin  is  dependent  upon 
the  net  heat  flux  conducted  into  the  fin.  The  boundary  condition  takes  the 
form: 


dT 


- K 


F 


x - 0 


q , = q I q - ci 

net  aero  ^part  rad 


(4) 


whe  re 


laero  = the  aerodynamic  heat  flux  at  the  fin  surface,  taking  into 
account  injection  of  the  fluid  into  an  air  boundary  layer 


‘Ipart  “ t^le  ^leat  flux  due  to  the  impinging  of  weather  particles 


clra(l  - the  heat  flux  due  to  radiation  from  the  fin  surface 

These  heat  flux  terms  may  be  evaluated  by  the  analysis  which  follows. 

Recause  some  percentage  of  the  coolant  vaporizes  at  the  fluid-air 
interface,  a boundary  layer  composed  of  air,  11,0  liquid,  and  11,0  vapor  is 


formed  over  the  TCNT  surface  as  shown  in  Fig.  1. 
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n 

R 

h 


Here  the  relationship  holds  that 


m = m + m , 

v 1 


whe  re 


m 


the  total  mass  flux  of  the  coolant 

the  mass  flux  of  the  coolant  which  vaporizes 


the  mass  flux  of  the  coolant  which  does  not  vaporize  at 
the  liquid-air  interface  and  therefore  enters  the 
boundary  layer  as  a liquid 


The  surface  temperatures  of  both  the  fluid  and  fin  are  sufficiently  low 
and  the  pressure  sufficiently  high  to  preclude  the  possibility  of  gas  phase 
reactions  at  the  surface.  Thus,  the  only  gas  phase  species  present  at  the 
surface  are  water  vapor,  molecular  nitrogen,  and  oxygen. 


The  width  of  the  fin  is  assumed  to  be  small  enough  such  that  the 
boundary  layer  does  not  relax  significantly  from  its  state  over  the  slot. 
That  is,  it  is  assumed  that  the  aerodynamic  heat  flux  to  the  fin  surface  is 
essentially  the  same  as  the  heat  flux  to  the  fluid  issuing  from  the  upstream 
slot.  The  only  difference  between  the  two  values  is  that  the  heat  flux  to  the 
fin  surface  is  corrected  to  reflect  the  higher  fin  surface  temperature.  In 
particular,  the  same  reduction  in  heat  transfer  due  to  vapo rization  at  the 
slot  surface  is  presumed  to  also  apply  at  the  fin  surface.  Thus,  the  aero- 
dynamic heat  flux  to  the  fin  surface  is  represented  by 


C 


q (fin  surface)  = q 
aero  no 


1 - 


H 


II 


R i 


C 


(5) 


H, 


c 

I i 


- 


i 


I: 


-14- 


whe  re 


qo  = the  cold-wall,  non-blown  aerodynamic  heat  flux  to  the 
surface  (evaluated  at  the  stagnation  point  in  this 
report! 

h = the  enthalpy  of  the  mixture  of  species  at  the  wall  surface 
(therefore  evaluated  at  the  fin  surface  temperature) 

H^  = the  recovery  enthalpy 

Cjj  = the  Stanton  number  for  heat  transfer  with  blowing 
(assumed  equal  to  the  Stanton  number  for  mass 
transfer) 


C 


ae  ro 


H 


0 u (H 
~e  e 


R 


h ) 

w 


C 


j,  = the  Stanton  number  in  the  limit  for  non-blown  heat 

o transfer  (m  =0) 
v 


C 


Ho  PeueHR 


The  terms  in  Eq.  (5)  may  be  determined  in  the  following  manner. 


The  ratio  of  Stanton  numbers,  CTT/CTI  , is  by  definition  a function  of 

H Ho  ’ 

the  mass  transfer  driving  force,  B1.  A commonly  used  expression,  with 

3 

some  theoretical  justification  , is  given  by 


'H  _ In  (AB ' -t  1) 


C 


H, 


AB' 


(A  = 1. 28) 


m 


B'  s 


p u c: 
e II 


(6) 


W.  H.  Dorrance,  Viscous  Hypersomc  Flow, 


McGraw-Hill  (1962). 
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The  ratio  of  Stanton  numbers  is  also  known  to  be  the  "blowing  correction" 

for  the  mass  flux  of  the  coolant  which  vaporizes  at  the  surface,  rh  , Follow- 

v 

ing  the  analysis  described  in  Appendix  B,  the  Knudsen - Lang muir  equation  for 
non-equilibrium  water  vaporization  is  incorporated  to  yield  the  following 
equation,  which  may  be  solved  iteratively  for  B'  for  given  values  of  fluid 
surface  temperature 


(Xb1  + i) 


T5  (p  ' ph2o)  = g/»ph2o 


The  variables  are  defined  in  Appendix  B. 

The  "wall  enthalpy"  h^  is  a function  of  the  enthalpies  of  the  individual 
species  present  at  the  surface,  so  that  we  may  write 


/_H_\  _»2°  /_H_\ 

VrtfL  28  VrtfA,  18  \rtf) 


where  the  species  mass  fractions  (K^)  existing  at  the  fin  surface  are  pre- 
sumed to  be  the  same  as  those  at  the  slot  surface  and  are  evaluated  in  terms 
of  the  mass  transfer  parameter  B'  (see  Appendix  B).  The  temperature  used  to 
evaluate  h^  is  that  of  the  fin  surface,  since  it  is  assumed  that  the  gas  phase 
adjacent  to  the  wall  adjusts  instantaneously  to  the  step  change  in  temperature 
from  the  fluid  value  to  the  fin  value.  It  is  recognized  that  this  instantaneous 
relaxation  in  temperature  is  inconsistent  with  the  previous  assumption  of  slow 
relaxation  with  respect  to  surface  blowing.  This  inconsistency  is  tolerated 
because  it  provides  a mechanism  for  the  fin  aerodynamic  heat  flux  to  be 
influenced  by  the  fin  surface  temperature  (which  seems  physically  plausible, 
at  least  to  some  degree).  A better  approach  would  be  to  use  weighted  averages 
of  surface  temperature  and  mass  transfer  parameter,  in  terms  of  relaxation 
distances,  to  determine  the  wall  temperature  and  blowing  corrections  to  the 


V 
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fin  heat  flux.  However,  in  light  of  the  boiling  heat  transfer  and  two-phase 
boundary  layer  uncertainties  in  the  analysis,  this  additional  complexity  does 
not  seem  warranted. 

The  parameters  HR  and  qQ  may  be  obtained  from  an  analysis  of  the 
boundary  layer  flowing  over  the  TCNT  nosetip.  In  this  report,  the  values  of 
recovery  enthalpy  and  cold  wall,  non -blowing  heat  flux  will  be  taken  as  those 
values  appropriate  to  a gaseous  air  boundary  layer.  That  is,  the  coolant 
mass  flux  entering  the  boundary  layer  as  liquid  water  m^  is  assumed  to  exert 
no  influence  upon  the  heat  flux  to  the  fin  surface.  This  is  obviously  incorrect, 
but  it  does  provide  an  upper  bound  for  the  surface  heat  flux  and  is  justified  on 
that  basis.  An  analysis  of  the  two-phase  air-water  boundary  layer  would  be 
required  to  accurately  define  the  heat  flux  and  is  beyond  the  scope  of  the 
present  effort. 

The  flux  at  the  fin  surface  due  to  particle  heating  q^  will  enter  the 
boundary  condition  [Eq.  (4)]  when  ambient  conditions  indicate  that  a "weather 
environment"  is  present.  For  particles  of  rain  or  snow,  the  particle  heat 
flux  is  a function  of  the  liquid  water  content  of  the  atmosphere 


01  0. 


(8) 


where 

Ci  = accommodation  coefficient,  which  is  experimentally 

determined  (herein  assumed  to  be  0.  7) 

P wr-  = liquid  water  content  with  units  of  mass  of  water  per 
volume  of  air 

U = free-stream  velocity 

oo 

Since  the  weather  test  case  that  shall  be  considered  later  does  not  include 
impingement  of  dust  particles,  Eq.  (8)  will  be  used  to  represent  particle 
heating. 
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The  radiative  term  in  boundary  condition  [Eq.  (4)]  usually  will  be 
much  smaller  than  the  aerodynamic  or  particle  heat  fluxes,  but  for  com- 
pleteness is  evaluated  as  follows: 


‘rad 


= £ <7T , 


(9) 


whe  re 

( = the  total  hemispherical  emissivitv  of  the  fin  assumed  to  be 

0.  16  (for  stainless  steel) 

a = the  Stefan- Boltzmann  constant 

Tp  = the  surface  temperature  of  the  fin 

C.  BACKFACE  BOUNDARY  CONDITIONS  FOR  THE 

FIN  AND  FLUID  EQUATIONS 

The  second  boundary  condition  used  for  the  fin  equation  and  the  boundary 
condition  for  the  fluid  equation  are  both  specified  at  the  backface  of  the 
respective  elements.  Both  elements  are  con-'idered  to  be  semi-infinite 
slabs;  that  is,  the  length  1 is  assumed  to  be  much  greater  than  the  thermal 
penetration  depth  of  either  fin  or  fluid  for  the  time  scale  of  interest.  Thus, 
the  backface  temperatures  remain  at  their  initial  values  for  the  entire  time 
considered,  so  that  our  boundary  conditions  are  formulated  to  be 


^F  ^F  initial  ^initial 


^B  ^B  initial  ^initial 


(10a) 

(10b) 


So  that  no  convective  heat  transfer  occurs  at  the  lower  portion  of  our  model, 
we  assume  that  the  initial  temperatures  of  the  fin  and  the  fluid  are  equal. 

Since  the  reentry  vehicle  is  in  storage  for  a considerable  length  of  time  before 
being  used,  the  fin  and  fluid  are  assumed  to  be  at  room  temperature  initially. 


! i 


i 


*s 


; 


-18- 


1 

J 


\n  alternate  formulation  of  these  boundary  conditions  is  to  consider  the 
conduction  in  both  the  fin  and  the  fluid  to  be  zero  at  the  backface.  This 
condition  is  perfectly  valid  for  the  semi -infinite  slab  assumption:  however, 
since  conduction  in  the  fluid  is  neglected  anyway,  it  provides  redundancy  as 
a backface  fluid  boundary  condition.  Since  we  wish  to  have  the  backface 
temperatures  for  both  elements  in  equilibrium,  the  relations  Eq.  (lOal  and 
Eq.  (lOb)  will  be  applied  to  the  problem. 

D.  GLOBAL  ENERGY  BA  LANCE 

In  order  to  check  the  validity  of  the  temperature  distributions  obtained 
by  solution  of  the  two  governing  equations,  an  overall  energy  balance  of  the 
system  should  be  satisfied.  We  consider  a control  volume  for  the  entire  sys- 
tem (composed  of  one  fin  element  and  one  fluid  element)  as  shown  in  Fig.  3. 


— a b »-| 


Fig.  3.  Control  Volume  for  Energy  Balance 
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The  net  energy  into  the  system  must  be  equal  to  the  change  in  sensible 
energy  in  the  fin  and  in  the  fluid,  yielding  the  relationship 


a / V:PF-3rLdx  + b J pbc 


l 


dT. 


f 


dT 


B 


'B^PB  dx  = a qnet  + 


+ b rnCpB  iTg 


dT. 


- T 

*-/  B 


x J - a kf"^T  x=£ 


(li) 


Equation  (11)  shall  be  used  to  test  the  accuracy  of  the  finite  difference 
solution  of  the  differential  equations. 
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III.  NUMERICAL  ANALYSIS 


ll 


The  two  coupled  partial  differential  equations  which  govern  the  behavior 
of  the  fin  and  fluid  temperatures  must  be  solved  by  numerical  methods.  An 
implicit  finite  difference  scheme  shall  be  used  to  represent  the  terms  in  the 
differential  equations.  An  explicit  formulation  was  initially  attempted,  but 
numerical  instabilities  caused  a quick  reversion  to  the  method  finally  em- 
ployed. 

The  basic  principle  underlying  an  implicit  finite  difference  representa- 
tion of  a partial  differential  equation  is  that  the  spatial  derivatives  are  repre- 
sented as  functions  of  the  dependent  variable  at  the  future  time.  Thus,  to 
obtain  spatially  dependent  values  of  the  variable  from  one  time  level  to  the 
next,  a system  of  equations  must  be  solved.  As  the  time  and  spatial  incre- 
ments used  in  the  difference  approximation  become  smaller,  the  values 
approach  the  true  solution. 

A.  FIN  EQUATION  AND  BOUNDARY  CONDITIONS  IN 

FINITE  DIFFERENCE  FORM 

In  order  to  ensure  stability  of  the  solution  even  with  severe  heat  fluxes 
at  the  fin  surface,  a variable  spatial  increment  will  be  incorporated  to  pro- 
duce a tighter  mesh  at  the  surfaces  of  both  fin  and  fluid.  The  three-point 
central  difference  representations  with  variable  mesh  size  are 

n + 1 R 2 Tn^ -1  + (1  - R2)Tn+1  - Tn+.1 

■ (1  + R)r4x. 

J J 


5T 

dx 


dzr 


n + 1 


RT 


n+1 
i + 1 


- (1 


+ R)Tn+1 
J 


+ t"^1 

ill 


Rm^x.)' 


whe  re 


j = the  spatial  node 


n = the  time  le ve  1 


R = A x.  ,/As,  the  ratio  of  successive  spatial  increments 
J-l  J K 

An  iterative  technique  shall  be  employed  to  determine  parameters  at  future 
time  levels.  Thus,  quasi- linearization  may  be  used  to  represent  more 
accurately  those  parameters  which  change  dramatically  with  time.  Quasi- 
linearized  formulations  of  these  parameters  are  derived  in  Appendix  C. 

The  resulting  implicit  finite  difference  representation  of  the  fin  equa- 
tion is  as  follows: 


pl  cl 

F PF  -t,  i + l ^ n ,,i 
~St  rF.  ■ TF.  = KF 


J+j 1 ls± 

Si  + r)R(4x.)2 

^ j 


/ aK{\  r2tfU1  +fI-R2,TFi  + 1 ' TFi  + 1 

\ FCPF  ds  J (1  + R)rJ3TT  1 ,+  aQj 


where  i denotes  the  level  of  iteration  of  the  set  of  equations  and 


°r  * f m 


n S"  - T- 

j j 


is 


i + l T i 

1.  B. 

J J 


which  is  derived  in  Appendix  C.  Iteration  of  the  system  of  equations  is  per- 
formed until  the  difference  between  successive  temperature  iterations  is  less 
than  some  established  convergence  criterion.  Since  some  of  the  properties  in 
Eq.  (12)  do  not  change  dramatically  as  a function  of  temperature  (specifically, 
the  properties  of  steel),  they  may  be  approximated  as  values  at  iteration 
level  i rather  than  at  i+l. 
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The  surface  boundary  condition  in  the  fin  may  also  be  represented  in 
finite  difference  form  because  of  the  conduction  term.  Yet,  if  we  assume 
that  the  first  node  (j  = 1 ) occurs  at  the  surface  of  the  fin,  the  temperature 
evaluation  at  the  node  j-1  has  no  meaning.  Thus,  a forward  difference  repre- 
sentation of  the  derivative  shall  be  used  for  this  boundary  condition  at  j = l: 


? i + 1 -,i  • . i 

R T„  - (1  + R)Tr  + (1  + 2R)T„l 
F3  F2  * 1 

(1  + R)Rdx, 


= q 


L + l 

net 


whe  re 


(13) 


• i + 1 
^net 


r<3q 


^ K+I-v, 


S^„et 


' i+1  „ i 

V B1  ' B1 


which  are  evaluated  in  Appendix  C. 

The  backface  boundary  condition  in  the  fin  is  represented  simply  as 


T 


i + 1 

Fj-J  MAX 


- T 


initial 


(14) 


whe  re  j = J MAX  re  fe  r s to  the  node  at  which  x = H . 

B.  FLUID  EQUATION  AND  BOUNDARY  CONDITION  IN 
FINITE  DIFFERENCE  FORM 

Since  the  two  equations  governing  fin  and  fluid  temperatures  are  coupled 
by  the  source  term  Q,  it  follows  that  the  time  and  spatial  increments  used  to 
approximate  the  fluid  equation  should  be  the  same  as  those  in  the  fin  equation. 
In  addition,  the  sudden  increase  in  Q (when  the  fin  temperature  exceeds  the 
fluid  saturation  temperature)  strengthens  the  coupling  further,  and  incites 
oscillations  in  the  temperature  distributions  with  time.  Hence,  quasi- 
linearization of  the  source  term  is  used  to  rectify  this  situation.  Spatial 
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oscillations  also  tend  to  occur  when  a central  difference  scheme  represents 
the  spatial  derivative  in  the  fluid  equation,  thus  necessitating  the  incorpora- 
tion of  a two-point  forward  difference  (i.e.,  upwind  differencing).  The 
resulting  implicit  finite  difference  representation  of  the  fluid  equation  is 


n1  c 1 

KB^PB  ^ i + 1 _ n 

' At  tb.  - tb. 

j j 


_ (mCPB  + PBCPB  s) 


T 1+1  - T 1 
B ...  B. 

J + 1 J 


^Qi  + 1 

b J 


whe  re  Q*  + ^ 
J 

parameters 


is  expanded  in  the  same  manner  as  for  the  fin  equation.  The 
Pfi  and  Cpg  are  evaluated  at  increment  level  i by  the  assumption 


that  fluid  properties  do  not  vary  greatly  with  temperature. 


While  the  use  of  the  upwind  difference  scheme  for  the  convective  term 
eliminates  the  spatial  oscillations,  it  is  only  first-order  accurate  spatially, 
whereas  the  central  difference  scheme  is  second-order  accurate.  However, 
sufficient  accuracy  may  be  ensured  by  using  a fine  mesh  near  the  surface 
(i.e.,  R < 1).  In  any  event,  excessive  concern  with  higher  order  accurate 
numerical  schemes  is  hardly  warranted  when  the  uncertainties  in  the  physical 
model  drive  the  solution. 


The  backface  boundary  condition  in  the  fluid  is  represented  simply  as 


T 


i+1 

Bj=JMAX 


T 


initia  1 


(16) 


C.  SOLUTION  OF  THE  SYSTEM  OF  EQUATIONS 

In  examining  the  system  of  Eqs.  (12)  and  (15),  we  find  that  the  two 
dependent  variables  are  evaluated  at  three  spatial  levels:  j-1,  j,  and  j + 1. 
Thus,  since  the  equations  can  be  written  in  tridiagonal  form,  the  unknown 
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future  temperatures  and  T_*'^  may  be  found  using  established  solution 

4 b j Bj 

procedures.  A system  of  tridiagonal  equations  takes  the  form 


A(j,  1)  u.+1  + A (j,  2)  m + A (j,  3)  u.  _ j + A(j,  4)  v.  + 1 + 

+ A(j,  5)  v.  + A(j,  6)  v j = A(j,  7)  (17a) 

and 

B (j,  1)  ujHl  + B(j,  2)  u.  + B(j,  3)  u + B ( j , 4 ) w+1  + 

+ B(j,  5)  v.  + B (j , 6 ) v._j  = B(j,  7)  (17b) 

where  u and  v are  dependent  variables  evaluated  at  time  level  i + 1,  the  j values 
represent  the  nodal  points,  in  this  case  varying  from  1 to  JMAX,  and  A and  B 
represent  coefficient  arrays  which  are  defined  in  Table  1. 

Clearly,  the  two  governing  difference  equations  can  be  represented  in 
this  manner  if  the  fin  equation  takes  the  form  of  Eq.  (17a)  with  T replacing  u 

r 

and  Tg  replacing  v and  if  the  fluid  equation  takes  the  form  Eq.  (17b).  The 
coefficients  in  the  equations  are  listed  on  the  following  pages. 

It  is  apparent  that  the  surface  boundary  condition  in  the  fin,  Eq.  (13), 
cannot  be  rearranged  into  tridiagonal  form  at  j = l since  it  contains  a variable 
at  the  j+2  level.  However,  if  solved  simultaneously  with  the  difference  equa- 
tion at  j = 2 , the  boundary  condition  can  be  manipulated  to  take  the  form  of 
Eq.  (17a).  The  backface  boundary  conditions  in  both  equations  are  easily 
represented  in  tridiagonal  form. 


4 

R.  D.  Richtmyer  and  K.  W.  Morton,  Difference  Methods  for  Initial  Value 
Problems,  p.  276,  Interscience  Publishers  (1967). 
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STABILITY  ANALYSIS 


To  ensure  that  the  difference  scheme  chosen  to  represent  the  equations 
will  yield  stable  solutions  for  fin  and  fluid  temperature  distributions,  a sys- 
tem stability  analysis  is  performed,  as  described  in  Appendix  D.  For 
various  phase  angles  unconditional  stability  requires  that  the  charac- 

teristic eigenvalues  (X)  of  the  difference  equations  are  less  than  or  equal  to 
one.  The  relations  for  the  two  eigenvalues  derived  in  Appendix  D have  been 
solved  as  a function  of  phase  angle,  assuming  a uniform  mesh  and  typical  fin 
and  fluid  properties  (i.e.,  constant  values  not  dependent  on  temperature). 

The  results  clearly  show  that  the  system  of  equations  is  unconditionally 
stable  for  all  phase  angles.  As  might  be  expected,  when  time  and  spatial 
increments  are  decreased,  the  eigenvalues  become  even  smaller.  The 
resulting  eigenvalue  for  various  time  and  spatial  increments  are  shown  as  a 
function  of  phase  angle  in  Appendix  D,  The  X versus  /?  curves  are  symmetric 
about  /3  = 180  deg. 


!■ 


s 


-27- 


■ — 


IV.  RESULTS 


Solutions  derived  from  the  analytical  formulation  are  obtained  and 
compared  with  experimental  data  in  this  section.  Two  cases  of  particular 
interest  are  considered.  First,  comparisons  are  made  with  surface  tem- 
perature data  obtained  from  a clear-air  ballistic  range  shot.  Following  this, 
calculations  are  presented  for  a typical  flight  test  trajectory  in  a weather 
envi  ronment. 

A.  BALLISTIC  RANGE  TEST  CASE 

The  test  data  for  the  ballistic  range  experiment  include  actual  TCNT 
surface  temperatures  located  in  the  stagnation  region  as  a function  of  time. 
Prescribed  coolant  mass  flow  rates  and  environmental  conditions  are  then 
supplied  to  the  numerical  solution  so  that  predicted  fin  temperatures  may  be 
compared  with  measured  stagnation  temperatures. 

The  data  shown  in  Fig.  4 indicates  that  the  surface  temperatures 
obtained  experimentally  are  far  greater  than  those  determined  numerically. 
Several  factors  could  have  caused  this  discrepancy.  The  most  predominant 
is  the  Thom  Correlation  which  represents  the  source  term  for  subcooled 
boiling  heat  transfer.  If  the  correlation  causes  too  much  heat  to  be  removed 
from  the  fin,  predicted  fin  temperatures  will  lie  below  the  actual  values. 
Another  case  was  run  to  establish  an  upper  bound  on  fin  temperatures  by 
using  a constant  Nusselt  number  for  the  source  term.  The  temperature  distri- 
bution for  Nu  = 4.0  (Fig.  4)  lies  above  the  experimental  data  points.  Thus, 
it  appears  that  because  of  the  influence  of  the  source  term,  a more  realistic 
correlation  for  the  heat  transfer  mechanism  should  yield  more  accurate 
results . 

The  possibility  exists  that  the  assumption  of  a laminar  boundary  layer 
also  causes  deviation  from  actual  surface  temperatures.  If  the  non-blown 
aerodynamic  heat  flux  is  greater  than  that  assumed  (i.  e.,  a turbulent  instead 
of  a laminar  heat  flux),  both  fin  and  fluid  temperatures  could  be  higher  than 
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Nu  = 4.0 


Fig.  4.  Fin  Surface  Temperature  for  AEDC  Ballistic  Range 
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predicted.  However,  when  a turbulent  non-blowing  flux  is  incorporated  into 
the  solution,  the  increase  in  fin  temperature  is  not  enough  to  reproduce  the 
observed  temperatures.  It  may  then  be  concluded  that  the  reported  discrep- 
ancies are  the  result  of  an  overly  optimistic  correlation  for  the  source  term. 

Even  if  a realistic  correlation  for  subcooled  boiling  were  incorporated, 
the  numerical  solution  for  fin  temperatures  would  be  higher  than  that  indicated 
by  ballistic  range  data.  The  difference  would  arise  from  the  assumption  of  a 
gaseous  boundary  layer.  It  may  be  postulated  that  the  liquid  particles  present 
in  the  boundary  layer  actually  reduce  driving  enthalpy  with  phase  change, 
lessening  the  effect  of  the  heat  flux  at  the  fin  surface  so  that  temperatures 
are  lower.  Had  this  type  of  result  been  obtained  in  the  numerical  solution, 
a type  of  accommodation  factor  would  have  to  be  included  to  represent  the 
aerodynamic  heat  flux  more  accurately. 

Since  coolant  mass  flow  rates  are  relatively  large  for  the  ballistic 
range  test  case,  the  convective  term  in  the  fluid  equation  may  become  much 
larger  than  the  transient  term.  A case  in  which  the  fluid  equation  is  assumed 
steady  state  has  been  run,  yielding  virtually  the  same  temperature  distribu- 
tions in  both  fin  and  fluid  as  those  for  the  transient  case. 

B.  FLIGHT  TRAJECTORY 

The  computer  code  was  applied  to  a conventional  Kwajalein  reentry 
into  a postulated  weather  environment  using  a representative  stagnation  point 
coolant  flow  rate.  The  mass  flux  of  the  coolant,  shown  as  a function  of  alti- 
tude in  Fig.  5,  exerts  a stronger  influence  on  the  temperatures  because  of 
the  relatively  large  range  of  values  it  assumes.  Fin  surface  temperatures 
are  seen  to  rise,  then  fall  off  as  altitude  decreases  as  shown  in  Fig.  6, 
similar  to  the  behavior  of  coolant  flux. 

Note  that  the  constant  Nusselt  number  case  would  predict  surface 
melting  in  Fig.  6,  presumably  setting  upper  and  lower  bounds  for  the 
actual  fin  temperature  distribution.  The  amount  of  predicted  erosion  in 
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Fig.  5.  Stagnation  Point  Ideal  Coolant  Requirements 
(Kwajalein  Mission) 
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the  nosetip  dramatically  contrasts  the  two  source  term  correlations,  as 
seen  in  Fig.  7.  When  the  fin  surface  temperatures  are  high,  as  in  the  case 
arising  from  a constant  Nusselt  number,  the  erosion  coefficient  is  greatly 
reduced,  thereby  increasing  the  erosion  rate  significantly.  Even  though 
the  amount  of  material  eroded  is  highly  unrealistic  for  an  actual  run  (nearly 
twice  the  length  of  the  original  fin  element  is  removed),  the  results  indicate 
the  strong  stability  of  the  solution  even  with  drastic  parameter  variations. 

The  results  also  suggest  the  need  to  increase  flow  beyond  the  given  schedule 
to  maintain  surface  temperature  below  the  material  melt  value  for  the 
constant  Nusselt  number  solution. 

Figure  8 provides  a comparison  of  fluid  surface  temperature  and 
saturation  temperature  as  a function  of  altitude.  The  incorporation  of  a 
non-equilibrium  vaporization  rate  to  evaluate  the  mass  transfer  driving  force 
B'  allows  fluid  temperatures  to  exceed  the  saturation  temperature.  While 
this  only  occurs  for  low  coolant  mass  flow  rates  at  high  altitudes,  it  is 
apparent  that  liquid  and  vapor  phases  can  be  present  at  the  surface  of  the  slot. 
The  assumption  of  only  single-phase  flow  in  the  analysis,  which  generates 
two  governing  equations  instead  of  three,  may  still  be  considered  plausible, 
in  light  of  the  relatively  short  time  in  which  vaporization  occurs. 

The  low  mass  flow  rates  for  the  coolant  prescribed  for  the  beginning  of 
the  flight  trajectory  could  indicate  a greater  importance  associated  with 
the  transient  term  in  the  fluid  equation.  However,  the  test  cases  exhibited  in 
Fig.  9 show  negligible  differences  between  fluid  temperatures  obtained  by  a 
steady-state  fluid  equation  and  those  found  by  transient  analysis.  Both 
ballistic  range  and  flight  test  runs  indicate  that  the  transient  term  in  the  fluid 
equation  may  be  justifiably  eliminated. 
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Fig.  9.  Comparison  of  Fluid  Surface  Temperatures, 
Transient  versus  Steady-State  Numerical 
Solution  for  Flight  Environment 
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V.  CONCLUSIONS 


The  thermodynamic  behavior  of  a discrete  injection  transpiration-cooled 
nosetip  has  been  modeled  by  considering  the  energy  exchange  between  the 
metal  matrix  and  coolant  for  a simple  fin  and  slot  combination.  The  relevant 
partial  differential  equations  governing  the  energy  transfer  in  the  fin  and 
coolant  have  been  derived  and  solved  numerically  using  an  unconditionally 
stable  implicit  finite  difference  scheme. 

While  the  heat  transfer  mechanism  between  fin  and  fluid  accounts  for 
the  possibility  of  subcooled  boiling,  it  is  maintained  that  the  bulk  of  the  coolant 
is  in  the  liquid  phase.  The  analysis  treats  the  governing  final  equation  as 
such,  but  in  order  not  to  preclude  the  possibility  of  two  phases,  a non- 
equilibrium vaporization  rate  of  the  fluid  into  the  boundary  layer  is  incorpo- 
rated. It  is  further  assumed  that  the  boundary  layer  of  the  nosetip  is  solely 
gaseous,  so  that  analysis  may  be  simplified  and  an  upper  bound  on  the  sur- 
face heat  flux  may  be  established. 

Comparison  of  nosetip  temperature  distributions  obtained  by  numerical 
means  and  by  experiment  indicates  that  the  correlation  for  subcooled  boiling 
heat  transfer  is  overly  optimistic  for  indepth  energy  exchange.  A more 
realistic  representation  might  yield  surface  temperatures  closer  to  the  actual 
values. 

While  both  governing  equations  are  considered  to  be  transient  throughout 
this  analysis,  a steady- state  postulate  in  the  coolant  equation  yields  virtually 
the  same  results  as  the  complete  transient  representation,  thereby  validating 
this  commonly  used  assumption. 
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APPENDIX  A.  DERIVATION  OF  GOVERNING  EQUATIONS 


[ 


A.  FIN  EQUATION 

The  physical  situation  depicted  in  Fig.  A-l  represents  a fin  volume 
element,  located  at  an  arbitrary  distance  z from  the  original  surface  of  the 
nosetip.  The  coordinate  x denotes  the  distance  of  the  element  from  the 
"actual"  surface  location  at  some  time  t,  that  is,  the  location  of  the  surface 
after  the  fin  has  begun  to  erode  at  a rate  s. 


Fig.  A-l.  Fin  Volume  Element 

The  T and  T-.  are  used  to  denote  the  temperature  of  the  fin  and  the  tempera- 
r 13 

ture  of  the  fluid  (or  "bulk"),  respectively.  Thus,  an  energy  balance  on  this 
volume  element  will  include  the  following  terms 

aTF 

% = KF  0x 


4 


| 


. 

; 


is  the  rate  of  heat  transfer  per  unit  area  due  to  conduction  in  the  fin  (where 
Kp  is  the  fin  thermal  conductivity). 


% = h(TF  - TB> 


is  the  rate  of  heat  transfer  per  unit  area  due  to  convection. 

The  change  in  sensible  energy  of  the  system  in  time  4t  is 

(PFCpF  • a . w • 4z)  • [tf  - tf  '1 

z , t |_  z , 1 4 At  a,  t 


so  that  an  overall  energy  balance  yields 
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- 2h  • w • 4z  • 4t  T - T 

F - - B - 
L z. t z, 


Dividing  by  a . w • 4z4t  and  taking  the  limit  as  4t  and  4z --►(),  we  obtain 
the  following  partial  differential  equation 


P 


F _ a 


stfV  .h 


fpf  at  az  v f az 


+ 2 a <TB  - V 


(A  - 1 ) 


In  order  to  incorporate  the  erosion  term  in  the  equation,  we  perform  a 
transformation  to  z-coordinates,  using  the  relation 


z = x + / s dt 


(or  dz  = dx  + s dt) 


(A-2) 


Our  goal  is  to  transform  temperature  T as  a function  of  z and  t to  T as  a 

F r 

function  of  spatial  coordinate  x and  time  t.  Thus,  if  we  wish  to  have 


TF(z,t)  = TF(x,t) 
we  may  then  differentiate  to  obtain 


*TF.  /^d,: 

~&rdz  +~dr dt 


dTr 

F 


dT, 


dx  + 


dt 


dt 


If  we  let  t = t and  dt  s dt,  then  using  Eq.  (A -2)  we  find  that 
6T  . dT  dT  dT 

~err  (dx  + s dt)  + dt  = dx  + dt 


or 


dt  dT 

F • . F 

~9T  s + dt 
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0T 
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dT, 


5x 


dx 


Equating  coefficients,  the  following  relations  result 

&Tf  _ 

dz  dx. 


dTF  dTF  . dTF 
dt'S 


Using  the  transient  heat  conduction  equation  for  the  fin  in  Eq.  (A-l),  we  obtain 


Pr  C 


0tf  d L 0tf'' 


f^pf  ~7TC  ~ dx  i f ~d 


+ P^c 


F PF 


dT 


F + 2 J (tb 
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(A -3) 


where  t now  replaces  t.  It  should  be  noted  that  when  no  erosion  takes  place, 
s equals  zero,  and  the  form  of  Eq.  (A-3)  reduces  to  that  of  Eq.  (A-l). 
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When  erosion  of  the  nosetip  surface  occurs,  examination  of  the  model 
indicates  that  the  entire  fin  of  length  / is  shifted  downward  by  the  amount 

rs  tit.  For  the  particular  application  of  the  transpiration-cooled  nosetip, 
the  length  £ is  defined  to  be  the  distance  from  the  surface  of  the  fin  (before 
erosion)  to  the  small  opening  through  which  coolant  is  injected,  as  shown  in 
Fig.  A-2. 


Fig.  A-2.  Control  Voluine  Geometry 

Thus,  if  the  fin  is  to  be  considered  to  have  constant  length  l after  erosion 
takes  place,  the  dimension  a in  Eq.  (A  -3)  must  be  replaced  by  a1  for  the 

lower  portion  of  the  fin  (i.e.,  for  x - / - / s dt  to  x =/ ).  The  same  back 

.° 

face  boundary  condition  would  apply  as  if  s were  zero,  since  we  are 
assuming  the  fin  to  be  semi-infinite. 


B.  FLUID  (BULK)  EQUATION 


A similar  energy  balance  may  be  performed  on  the  fluid  volume  ele- 
ment shown  in  Fig.  A-3,  with  x-  and  z-coordinates  identical  to  those  of  the 
fin  element  so  that  erosion  may  be  taken  into  account. 


-I , \ I | FIN  AND  FLUID  SURFACE 
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Fig.  A-3.  Fluid  Volume  Element 


3t 


B 


The  terms  in  the  energy  balance  are 

= 'kb 

is  the  rate  of  heat  transfer  per  unit  area  due  to  conduction  in  the  fluid. 


qc  = MTf  - Tb) 


is  the  rate  of  heat  transfer  per  unit  area  due  to  convection  over  the  fin. 


mCPTB 


is  the  rate  of  heat  transfer  per  unit  area  due  to  injection  of  the  fluid,  where 
m is  the  specified  mass  flow  rate  of  the  coolant  per  unit  area. 
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An  overall  energy  balance  yields  the  relation 
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from  which  the  partial  differential  equation  is  obtained 


15  PB  7TT 


‘B  dFj 


+ ('hcPB»  ST  + 2 ■ S (Ti 


(A  -4) 


A coordinate  transformation  similar  to  that  for  the  fin  equation  yields  our 


final  result 


n n 

pi}cpB~&r 


= ^ (KB  ) + (PBCPBS  + m(  Pb)  751 T + 2 b (TF  ' rB' 


(A  - 5) 


/Because  erosion  causes  the  fin  to  be  "shifted"  downward  by  the  amount 

s dt,  the  region  of  interest  in  the  fluid,  also  of  length  /,  must  translate 

downward  as  well.  Referring  to  Fig.  A-2,  it  is  apparent  that  the  dimension  b 

in  Eq.  (A-5)  must  be  replaced  by  b1  for  x - l - I s dt  to  x s /.  It  is  clear 

that,  since  m is  defined  as  the  mass  flow  rate  of  the  coolant  per  unit  area,  a 

different  value  rn1  should  be  substituted  into  the  equation  for  / - / s dt  < x < f . 

• • **  o 

From  continuity,  this  new  mass  flow  is  found  to  be  m'  = m(b/b'). 

In  the  derivation  of  Eq.  (A-5),  the  fluid  was  assumed  to  be  incompres- 
sible such  that  work  due  to  volumetric  changes  is  negligible.  The  fluid 
velocity  is  assumed  constant  such  that  viscous  dissipation  can  be  ignored. 
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APPENDIX  B.  VAPORIZATION  FROM  A LIQUID  WATER  SURFACE 


1 


The  mass  transfer  parameter  B1  is  required  to  determine  the  reduction 
in  aerodynamic  heat  flux  due  to  coolant  vaporization  at  the  slot  surface.  The 
mass  transfer  parameter  is  obtained  from  the  analysis  given  in  this  appendix. 


The  problem  is  to  determine  the  vaporization  from  a Liquid  water  sur- 
face over  which  an  air  boundary  layer  is  flowing.  In  order  to  simplify  the 
analysis,  the  usual  assumptions  of  unity  Prandtl  and  Lewis  number  and  equal 
diffusion  coefficients  are  invoked.  With  these  assumptions,  the  elemental 
mass  fractions  at  the  surface  are  obtained  from  the  wall  boundary  conditions 
for  the  elemental  conservation  equations  in  terms  of  the  mass  transfer  para- 
meter. 
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(B-l) 


where  subscripts  w and  e refer  to  wall  and  boundary  layer  edge  conditions, 
respectively,  and  superscript  tilde  refers  to  element  mass  fractions. 
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Assuming  that  the  surface  temperature  is  low  enough  to  preclude 
dissociation  reactions,  the  gas  phase  at  the  wall  consists  solely  of  molecular 
oxygen,  nitrogen,  and  steam.  Elemental  conservation  at  the  wall  yields 


K +i^K 
Oz  + 18  ^HzO 


KH  = 18  KH,0 
w 2 


(B-2) 


The  sum  of  the  partial  pressures  must  equal  the  mixture  pressure 


PNZ  + po2  + pH2o  = p 


The  partial  pressures  are  related  to  the  species  mass  fractions  as 
follows : 


(B-  3) 
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(B-4) 
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If  it  is  assumed  that  the  vaporization  reaction  is  in  equilibrium  such  that 
the  water  vapor  partial  pressure  is  identical  to  its  equilibrium  vapor  pres- 
sure (i.e.,  Ppj  q = P^j  q),  the  mass  transfer  parameter  maybe  obtained  from 
the  foregoing  expression  as 


,8^pH2o 


p - p; 


kn  ko 

o _ e , e 

p ~ 28  32 


(B-5) 


where  was  defined  for  notational  convenience. 


Since  the  vapor  pressure  is  a known  function  of  temperature,  the 
equilibrium  value  of  B*  is  easily  computed  for  specified  values  of  surface 
temperature  and  pressure.  The  equilibrium  mass  transfer  parameter  is 
undefined  for  surface  temperatures  greater  than  the  saturation  temperature 
(i.  e.,  P^,  „ > P).  Thus,  the  use  of  the  equilibrium  assumption  would 

h2o 

constrain  the  solution  for  the  surface  coolant  temperature  to  always  be  less 
than  the  saturation  value.  In  some  instances  this  is  physically  unrealistic, 
and  the  mass  transfer  parameter  will  be  obtained  from  a non-equilibrium 
a s sumption. 

The  non-equilibrium  vaporization  rate  is  modeled  using  the  Knudsen- 
Langmuir  equation,  which  is  given  as 


m = a 


p - p 
H2°  H2° 


(B-6) 


where  Ct  is  the  vaporization  coefficient,  and  M is  the  molecular  weight  of 
water  vapor. 

Using  the  definition  of  B'  and  dividing  through  by  the  heat  flux  yields 


G PH2C-PH,Oy 


P U Cu 
*e  e H 


(B-7) 


This  equation,  along  with  the  previously  obtained  expressions,  yields  a 
nonlinear  algebraic  equation  for  B1. 


p - pH2o ) = g*ph 


(B-8) 
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Since  the  ratio  of  blowing  to  non-blowing  Stanton  number  is  presumed  to 

be  a known  function  of  B',  this  equation  may  be  solved  iteratively  for  B1  for 

specified  values  of  surface  temperature,  pressure,  and  non-blowing  heat 

flux  (i.  e.,  peu eCjj  ).  The  Stanton  number  expression  used  in  this  report  is 
o 


given  as 


CH  _ /n(l  + XB') 
CH  XB' 

o 


A = 1.28 


As  equilibrium  is  approached  (i.e.,  Cj^  — >-0),  Eq.  (B-8)  yields  exactly 
the  value  of  Eq.  (B-5)  for  B1.  Furthermore,  Eq.  (B-8)  is  well  defined  at 
the  saturation  temperature,  thereby  allowing  solutions  for  surface  coolant 
temperature  in  excess  of  the  saturation  value. 

5 

The  required  value  of  vaporization  coefficient  was  obtained  and  is 
given  as  01  = 0.  04. 

Some  typical  values  for  B1  are  shown  in  Fig.  B-l  for  a pressure  of 
100  atm  for  both  the  equilibrium  and  non-equilibrium  case. 


John  P.  Hirth, "Evaporation  and  Sublimation  Mechanisms,  " The  Characteri- 
zation of  High  Temperature  Vapors,  Chapter  15,  ed.  T.  L.  Margrave, 

John  Wiley  and  Sons,  Inc.,  New  York  (1967). 
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APPENDIX  C.  QUASI  LINEARIZED  FORMULATION  OF  PARAMETERS 


The  two  parameters  in  the  governing  partial  differential  equations 
which  are  strong  functions  of  temperature  are  Q,  the  coupling  heat  transfer 
source  term  in  both  equations,  and  q , the  fin  surface  heat  flux  which 
enters  into  one  boundary  condition.  When  the  governing  equations  are  repre- 
sented in  finite  difference  form,  it  is  necessary  that  these  two  terms  in 
particular  are  evaluated  at  "future"  time  levels  to  ensure  consistency.  Hence, 
a technique  is  incorporated  to  evaluate  the  terms  iteratively.  The  terms  Q 
and  q will  be  represented  at  future  time  levels  by  the  method  of  quasi - 
linearization,  which  consists  of  simply  expanding  the  terms  in  a Taylor 
series. 

A.  THE  SOURCE  TERM  Q 


The  source  term  evaluated  at  iteration  level  i + 1 is 


: 


Since  Q is  represented  by  two  different  relations,  dependent  upon  the  regime 
of  the  fin  temperature,  the  derivatives  take  two  different  forms: 

IfTrisTsat(fluid» 
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where  P is  in  psia; 


= 2(0-0536)e(P/630)(TFj  ' Tsat)  ’ 


and 


(C-3) 


When  the  fin  temperature  exceeds  the  fluid  saturation  temperature,  the  rela- 
tion which  yields  the  maximum  value  of  Q is  chosen  to  represent  the  source 
term.  In  this  case,  the  associated  derivatives  for  the  particular  Q relation 
are  substituted  into  Eq.  (C-l). 


B.  THE  SURFACE  HEAT  FLUX  qnet 


The  net  heat  flux  impinging  on  the  surface  of  the  fin  may  be  divided  into 
three  major  components.  These  are  (a)  the  fluxes  due  to  aerodynamic  heating, 
(b ) particle  heating,  and  (c)  radiation.  The  representation  of  the  surface 
boundary  condition  is 
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The  particle  heat  flux  is  independent  of  either  temperature,  but  both  q and 

aero 

q ra(j  enter  into  the  evaluation  of  the  derivatives 
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To  obtain  the  derivatives  of  aerodynamic  heat  flux,  it  must  first  be 
determined  how  the  term  is  dependent  on  temperature.  For  the  relation 


q = q 
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the  terms  qQ  and  are  functions  of  the  local  boundary  layer  conditions  and 
therefore  independent  of  T or  T„,  the  term  h is  found  by 


It  is  apparent  that  h^  is  a function  of  the  fin  temperature  T^,.  Since  the 
mass  transfer  drivi  g force  B'  is  solved  iteratively  using  the  relation 


whe  re 


12 8.  0536a 


It  can  be  seen  that  both  h and  Cu/C  are  functions  of  the  fluid  temoerature 

W H Hq 

Tg.  Thus,  the  aerodynamic  heat  flux  derivatives  become 
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Ail  of  the  derivatives  in  Eqs.  (C-4)  and  (C-5)  may  be  evaluated  in  closed  forn 


APPENDIX  I).  STABILITY  ANALYSIS  FOR  SYSTEM  OF  EQUATIONS 


A.  GENERAL 


A stability  analysis  for  the  coupled  systems  of  difference  equations 
obtained  in  Section  III  will  provide  information  about  the  constraints  (if  any) 
imposed  on  selection  of  time  and  spatial  increments.  For  a system  of 
equations  to  be  unconditionally  stable,  it  is  required  that  the  characteristic 
eigenvalues  (X)  be  less  than  or  equal  to  unity  for  all  values  of  the  coefficients. 
I he  analysis  herein  has  assumed  the  coefficients  in  the  governing  equations 
to  be  extreme  "worst  case"  values,  and  the  spatial  increments  to  be  uniform 
(R  = 1). 


The  fluid  and  fin  difference  equations  may  be  represented  as  follows  for 
stability  considerations. 


1. 


FLUID  EQUATION 
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(D-l) 
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2. 


FIN  EQUATION 


T n-ri  _ T n _ c T n+1  . 2T  n+1  + T n+ 1 
j j J+l  J j-1 
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Since  constant  coefficients  in  the  difference  equations  are  assumed,  the  need 
for  iterations  is  eliminated,  and  the  unknown  temperatures  are  to  be  found  at 
time  level  n + 1. 

4 

In  accordance  with  Richtmyer  and  Morton  , it  may  be  assumed  that  the 
set  of  solutions  to  these  simultaneous  equations  takes  the  form 


n = Gn€J 

J.  U1 


g!?  € J 


where  G™  and  G™  are  some  functions  of  x and  t raised  to  the  n^  power  and 
£ = cos^  + i sin  ft,  a function  of  the  system  phase  angle  ft  (£  is  raised  to  the 
power  j). 
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Plugging  these  solutions  for  temperature  into  Eqs.  (D-l)  and  (D-2),  and 
using  the  identities 


and 


f - ^ = 2i  sin 
€ + j = 2 cos  /3 


the  following  is  derived  for  the  fluid  and  fin  equations 

Gj+^[l+B-A(cos/?  + i sin/3  - 1 )]  - G^ 

°2  = B 

n + 1 G2  + 1 t1  + E ' 2C(cos/3  - 1)  - i(2D  sin/3)]  - g£ 

G1  = E 


(D-3) 


(D-4) 


These  equations  may  be  solved  simultaneously  to  obtain  relations  for  the 
functions  G^  and  G^  at  time  level  n+1  in  terms  of  the  functions  at  time  level  n 
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with 


R s [1  + B - A (cos/?  - 1)]  [1  + E - 2C(cos/3  - 1)]  - (A  sin/?)(2D  sin f?)  - BE 

I s (A  sin£)l  1 + E - 2C (cos - 1)]  + [1  + B - A(cos/3  - 1)]  (2D  sinfi) 

The  matrix  containing  the  a^  coefficients  is  known  as  the  amplification 
matrix.  It  is  the  eigenvalues  (A)  of  this  matrix  which  yield  the  important 
information  concerning  system  stabilization 


all  ‘ * a 12 


21 
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The  final  polynomial  to  be  solved  for  the  eigenvalues  is 
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(D-5) 


whe  re 

xx  = R |2  + B + E - (2C  + A )(cos/?  - 1)  | - I(2D  + A)sin/S 

and 

x ^ = -R (2D  + A)sin£  + I |2  + B + E - (2C  + A)(cos0  - 1)| 

If  | A I ^ 1 for  all  values  of  A,  B,  C,  D,  and  E,  and  for  0°</J  < 360°,  then 
it  may  be  concluded  that  the  system  is  unconditionally  stable.  Since  extreme 
constant  parameters  in  this  system  of  equations  has  been  assumed,  the  terms 
A - E may  be  varied  by  changing  the  time  and/or  spatial  increments. 
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For  the  case  where  fl  = 0 , Eq.  (D-5)  reduces  considerably,  and  yields 
the  two  eigenvalues 


\ — 1 

1.2  “ 1 + B + E ’ 


Since  the  coefficients  B and  E are  always  greater  than  zero,  the  stability 
criterion  is  unconditiona lly  satisfied  at  fi  = 0°. 

When  * 0°,  the  solutions  for  A are  found  by  complex  analysis.  The 
resulting  eigenvalues  are  then 


^l  +Xz)  , 2 s'r 

(r2  + I2)2  R2  + I 

(x"  + X2  J - r- 

1 _1 LA  + r + 2 >,r. 

2 2 2 2 

R + I R + ] 


z \<r  / 

1 2 (X1 

R^  + I V 1 


6 , .6 

cos  ^ 


^ (X1 


0 , 0 

cos  2 + x2sin  j 


(D-6) 


whe  re 


r = Aa2  + b2 


a f - 1 b 
0 = tan  — 

a 


2 2 

X1  " X2  4 R 

Jr2  + i2)2  ' r2  + i2 


2x1x2 


(r2  + I2 


)2  ' K2 


ancl  Xj,  x,,  R,  and  I are  the  same  as  previously  given. 


When  combinations  of  At  and  Ax  are  substituted  into  these  relations 


for  /?  varying  between  0 deg  and  360  deg,  a complete  picture  of  the  system 
stability  is  obtained.  In  every  instance,  the  computer  solution  shows  that 
each  eigenvalue  is  less  than  unity,  except  at  fi  - 0 , when  becomes  equal 
to  one.  The  case  when  both  eigenvalues  approach  unity  closest  is  for  a very 
small  time  increment  and  a very  large  spatial  increment  (dt  = 0.00001  sec, 

4x  = 0.  1 ft).  Yet  the  eigenvalues  always  remain  less  than  or  equal  to  unity, 
indicating  that  the  governing  difference  equations  are  unconditionally  stable. 

A sample  plot  of  A versus  0is  shown  in  Fig.  D-l. 
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NOMENCLATURE 


B1 

b 


H 


R 


(H/RT). 

K 

K. 

1 

m 

R 

s 

T 

t 

l 

P 

P u cu 

e e H 

P u c 

e e Ho 

Subscripts 

B 

F 

W 


fin  width 

mass  transfer  parameter  (defined  in  Eq.  (6)] 

s lot  width 

specific  heat 

recovery  enthalpy 

normalized  enthalpy  of  species  i 

thermal  conductivity 

mass  fraction  of  species  i 

coolant  mass  flow  rate  per  unit  area 

universal  gas  constant 

fin  surface  recession  rate 

temperature 

time 

s lot  depth 
density 

cold  wall,  blowing  heat  flux  parameter  = q/H 

K 

cold  wall,  nonblowing  heat  flux  parameter  = q^/H^ 


coolant  properties 
fin  properties 
wall  values 
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